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1. Introduction 



It is well known that multipath fading can significantly effect the link reliability in a 
communications system or target detectability in the case of radar. The path between 
a transmitted and a receiver is often obstructed by natrual or man-made obstacles 
such as hills, buildings, atmospheric layers, trees, rain, fog, etc. Propagation outage 
due to multipath fading depends in a complicated manner on propagation climate, 
terrain features, path length, radio frequency, and fade margin. In the case of atmo- 
spheric multipath fading, interference due to two or more super- refracted rays arriving 
at the receiver via different paths can lead to a complete loss of signal. Other patho- 
logical phenomenon such as obstruction fading (caused by sub-refractive atmospheric 
effects) and ducting (caused by extreme super- refractive effects) are also possible. 
Such phenomenon are more common in warm and tropical climates, particularly near 
shores, where elevated inversions are formed easily due to the large temperature and 
partial pressure differentials. Reflection multipath fading, which is due to interfer- 
ence between the direct and the ground reflected ray depends strongly on the terrain 
geometry and ground constants. Moreover, elevated terrain features could completely 
mask a receiver from a transmitter leading to severe loss of signal (in some cases, it 
is advantageous to site antennas behind hills to provide shielding against undesirable 
interference). It is very important to assess the effects of environment on the link. A 
computer model that can take into account a given refractive index profile, terrain 
elevation data, and varying ground parameters will be very helpful in predicting the 
link performance. 

In this report we present formulation details for an efficient numerical solution of 
wave propagation in an inhomogeneous atmosphere and over irregular terrain using 
parabolic equation. 

Unlike all other previous formulations of the parabolic equation, we will use a modified 
Helmholz equation for propagation in an inhomogeneous atmosphere as suggested by 
Maxwell’s equations (all previous formulations use a Helmholtz equation which is 
only true for fields in a homogeneous medium). Because the parabolic equation is a 
full-wave method, it will include all aspects of wave propagation such as reflection, 
refraction, diffraction, and surface wave propagation. In this respect it is far superior 
to the commonly used ray method. 

Parabolic equation approximation to an elliptic partial differential equation, which 
the true fields satisfy, has proven to be a viable approach for studying propagation 
problems in underwater acoustics. The method is just gaining popularity with the 
electromagnetic community. Although the parabolic equation regards waves as es- 
sentially traveling one-way, it allows a rapid solution of the fields by way of marching 
along the range starting from an initial range. Another advantage of the PE method 
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compared to the ray methods is that it is valid even in the shadow region where the 
simple ray methods completely break down. Furthermore, it appears to be the only 
practical method for predicting propagation over long ranges (greater than 1 km) over 
a wideband from HF (a few MHz) through SHF (a few tens of GHz). The method is. 
however, not without limitations. In its standard form, the accuracy of the method 
is limited to waves traveling essentially within ±10° from horizontal. Furthermore, 
treatment of the boundary conditions on the uneven terrain is difficult. 

The method we propose to use will attempt to remove both of these defficiencies. 
Firstly, we use a Helmholtz-like elliptic equation to describe the fields and arrive 
at a wide-angle parabolic equation subject to certain approximations. To facilitate 
imposition of the boundary conditions on the irregular terrain, the equations will be 
transformed to a body-fitted curvilinear coordinate system. The PE will be solved 
using finite differences on a non-rcctangnlar mesh. This is a major departure from 
previous apj)roaches which will not only make the method more efficient but also 
more accurate. 

In Section 2, we derive the exact equations satisfied by 2D fields in an inhomogeneous 
atmosphere. Impedance boundary conditions are used to characterize the ground. In 
Section 3, we present details on the impedance boundary conditions. Starting from 
the exact equations presented in Section 2 for the fields, we derive, in Section 4, a 
parabolic equation (PE) valid for narrow angle propagation. This case is termed as 
the standard PE. The standard PE is generally valid for propagation angles that 
are within ±10° from horizontal. To accomodate waves traveling at larger angles, we 
present the derivation of a wide angle PE in Section 5. To truncate the computational 
domain we derive boundary conditions on an upper boundary, which are termed as 
the tropospheric boundary conditions. The derivation is based on the solution of 
Schrodinger type parabolic equation in a quarter plane t > 0. y > 0. This is presented 
in Section 6. 

Eor an efficient numerical implementation, we transform the differential equation and 
boundary conditions to a curvilinear coordinate system. This is presented in Section 
7. Details on the numerical generation of the curvilinear coordinate system are given 
in Section 8. Finally, in Section 9, we present steps leading to a finite difference 
solution of the equations using a Crank-Nicholson implicit scheme. 
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2. Solution of 2D Fields in an Inhomogeneous Medium 



Consider an electric source producing fields in an inhomogeneous region as shown 
in the figure below. Let us assume that both the sources and the medium are two- 
dimensiona) in nature in that all quantities are independent of the r-coordinate. As 
in the case of a homogeneous medium the fields can be decomposed into a case 
(vertical polarization) and a TM^ case (horizontal polarization). It is eissumed that 
propagation takes place in the xy-plane. 




V X F = (1) 

V X // = jijUcE 4- J (2) 

In the case of vertical polarization the, fields could be written in terms of the z- 
component of the magnetic field, 7/ = zH^ {TE^ fields). Substituting into (2) we 
have 

y//, X z - / 
juj 



-jujfiH 



V X // = V X (z//,) = VIE X z = jutE + J 
In a source-free environment, we have 



tE = 



eE = 



VH, X z 



Substituting into (1) we get 

_ - 1 _ /V//, z _ /V//, \ 

V X E = T-V X X z = -—V • = 

ju> \ e J JUJ \ € / 
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The equation satisfied by the magnetic field is then 






V- (^) = 0 



Vertical Polarization 



Once the magnetic field is determined the electric field is given by 






i X V//, 



Note that Hz does not satisfy the Helmholtz equation unless t is constant. 
For horizontal jjolarization on a similar analysis with E = zE^ shows that 




Horizontal Polarization 



(3) 

(4) 



(3) 



If the medium is non-magnetic, /i = piQ and satisfies the Helmholtz equation. W’e 
will characterize the ground in terms of impedance boundary conditions [3]. 



\\e may combine the vertical and horizontal cases shown in (3) and (5) into an 
equation of the form 

V ■ (qVv) + Bv - 0 (6) 

where 



Q 



V; 




TE or \'ertical Pol. 



— TM or Horizontal Pol. 

akln^{x,y) > 0 

[ H, TE Pol. 

[ E, TM Pol. 



(U 

( 8 ) 

(9) 



The quantity n{x,y) is a position dependent refractive index of the medium. The 
partial differential equation (PDE) given in (6) is elliptic, for we have 



\ dx J 



dy (“aj) 



-f /3xl> = 0 
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or 



Q 






dv da dll- da 

+ — T — ~a — ~ 0 

dx dx dy dy 



dx"^ dy^ 

The discriminant for the above PDE is —1 <0 implying elliptic nature of the equation. 
Equation (6) could be expressed as 



1 da 



1 da 



VV’ + - — vv + - + -V = 0 

a ax a ay 



0 , 
— i 

a 



Eor a non-magnetic. loss-less medium, we have 

1 



a — < 



CoCr 



Vertical Pol. 



— , Horizontal Pol. 
^>0 



so that for vertical polarization 



1 da 
a dx 






dx \Cr 
2 ^ 
n dx 



1 dn‘ 



dx 



{n is the refractive index) 



Similarlv. 



Letting 



1 5q 2 dn 

a dy n dy 



Qii^-.y) - 



a2{x,y) = < 



2 dn 
n dx 

0 

2 ^ 

n dy 

0 



\^ertical Pol. 
Horizontal Pol. 
Vertical Pol. 

Horizontal Pol. 



equation (10) may be expressed in the form 

VV + + a2‘>Py + kln'^tp = 0 



( 10 ) 



( 11 ) 
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3. Impedance Boundary Condition 



Impedance boundary condition relates the tangential components of electric and 
magnetic fields at the interface of two media. If r/ is a unit normal and s is a unit 
tangent as shown in the figure below, the boundary conditions are given by [3] 




s = X cos 6 y smO 
V = —X sin 0 y cos 0 

X = s cos 6 — 1/ sin 9 
y = s cos 6 + 6 cos 6 

( 12 ) 



where = ZJrjo is the surface impedance normalized to the free space value t]q. 
The equation may also be expressed as 



(;/ • E )6 — E = —t)oAsI/ X // 



1 / X 



E = r/oAjf- X (i> X fj) 



= 7?oA, {v ■ H )6 - H 



( 13 ) 



The surface impedance is determined from the intrinsic impedance of the medium 
by considering plane wave reflections from the interface. The complex propagation 
constants, 71, 72, and the intrinsic impedances 7/1 and 72 in terms of the media 
constants arc indicated in the figure. 



'll — i<^V^riPo(o'i + JU'eoCrl) 




= - J—'] 

^ol^T\^rc\ 

72 = -^q6t2(-tc2 




According to Snell’s law, we have 71 sin 6 ^ = 72 sin 6 t 



The plane wave reflection coefficients for the vertical and horizontal polarizations, Ry 
and Rm are [7] 



Rh 

Rv 



7/2 sec 6 t -Tji sec e , ^ ^ r t j ' tH 

burlace Impedance Z = 7/2 secy, 

7/2 sec Ot + Vi sec u, 

7/2 cos 0t - 7/1 cos 0, y 

a T =02 cos Oi 

7/2 cos y, — 7/1 cos y, 



A" = 



02 sec y, _ 1j2 
01 01 



1 — f — sin y, 



-1/2 



I^lr2 ^rcl 


1 — cos Ipi 


^ ^rc2 


/^r2^rc2 



and 



A = 



1 - 1/2 



/ /^t- 2 ^rcl 
flri ^rc 2 



1 COS 



I^ t 2 ^ tc 2 



1/2 



For the special case of /7i = //2 = yo, (T’i = 0, 



A = 



For normal incidence v, = 90° 





/ ^rl 


1 - 


67-1 


COS^ ll>i 




V ^r2 -J(^t 2 


^r2 ~ jTrr2 




1 


1 - 


^rl 


COS^ 7/), 




V Cr2 -j<^r 2 


^r 2 ^ j^r 2 



-1/2 



1/2 



(14) 

(15) 



A"^ = A'' = 



^rl 



^r2 “ J<^r2 

For the 2-D case, impedance boundary conditions for vertical polarization can be 
simplified as 

0xE = 7/oA]' [(i> • H)u - fl] = -t/oA^// 

Taking a dot product on both sides with z 
z ■ {u X E) = -tjqAI'IE 

Since z X {> ■= — s , we have 
-s ■ E = -t/oA^/Zj 



Substituting from equation (4) for E, we get 

. zxVH 



ju'e 



= -t/oA://, 
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0 ■ V//^ = jucTfoA] = jkotrA] 
i.e., the above can be put conveniently in the form 

IBC \^ertical Pol. (16) 

A similar analysis for horizontal polarization yields 

IBC Horizontal Pol. (17) 

This may also be obtained by resorting to duality. 

For a perfectly conducting material, we have = 0 and 

o r r 

- -- = 0 \ertical Pol. 

uu 

= 0 Horizontal Pol. 



-i ^ r - n 



dIE 

di' 



- jkotrA'^ IE = 0 
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4. Standard PE Derivation 



We will make some approximations and cast (11) in the form of a parabolic equa- 
tion which permits a rapid numerical solution. 






Let 'ipix.y) = — 






Then 



Vx = Ux - jkoV - 



u \ e 






2x^/2 / ^/x 



i^x -jkoy)- 



^-jkox 






X — ^ oc 



-jkox 



= ^y- 



\fx 



* ^yx — ^xy (^xy J^O^y)' 



7T 



^vv = — ■ ^xx ~ {^xx - ‘^jkoUr - klu)- 



s/x 



'yy - -yy ^ 

Substituting into (11) we get 

Uxx + J^yj, - ‘^jkoUj + OiUj. - 2jkoaiU + a2Uy + (n^ - l)^'oU = 0 

or 

Hxx + ^yy + (Qi ~ 2jT.'o)Ux + ^2^^y + — \ — 2] — ) A'qU = 0 

If now we impose the approximation that 

,1/2 



ji^xxl T 4 /tq^ |i/j 



we obtain 



-1 



Ur = 



(uj — 2jiA’o) 



Uyy + a^Uy T ( H ^ ~ 1 ~ 2 j ^) u 



or 




(18) 



This is the exact form of narrow angle PE approximation. W'e would also like to 
express the impedance boundary condition in terms of the ‘i/’ functions. 



X = SCOS 0 — {)smO 

dx ^ ^ ^ 

Xly = — = ly * Vx = U ' X 

dv 

= — sin0 

V = —X sin 0 + y cos 6 
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F^or vertical polarization we have 



dih 

du 



-iV.Ayy/, = 0 



dH, 

du 



- jkon^A''^H, = 0 



lh{x.y) 

dll, 

du 



-r^u[x.y) 

y/x 



- jkox^u 



ux^ \ 

2x3/2/ y/^ 



(Ui/ 



jkox^u) 



^-}koi 



X — ► cc 



• ■ - jkQ{n‘^A] + x^)u = 0 

Similarly for horizontal polarization, we have 



IBC \’ert. Pol. 



Uu - jko 




IBC Horz. Pol. 



We combine the two by defining 



' —jko{n^A]' — sinO) Vert. 

Horz. 



and writing as 



+ C]U = 0 



(19) 



The parabolic equation given in (18) is valid for propagating angles close to horizon- 
tal (±10° in practice) [l]. To accomodate waves at higher angles we would need a 
wide angle parabolic equation whose derivation is accomplished through a pseudo- 
differential operator formalism [1]. 
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5. Wide Angle PE Derivation 



Let us assume that the media constants are independent of horizontal range so that 
a{x,y) = q(j/). n{x,y) = n{y). We then have from (17) 



Uxx + ^yy ~ H ~i~^y + = 0 

a dy 



Let 



P = ^ , Q = 

ox 






1 52 



+ 



\ da d 



+ n‘ 



kl k^a dy dy 

Using this notation, the above PDE may be expressed as 



Pseudo differential 
operator in y 



[P^ -2jkoP + {Q^ -l)kl\u=0 

which we may factorize as 

{P - jko - jkoQ){P - jko + jkoQ)u = 0 
To see this we expand the operator on the left hand side to get 

P^ - jkoP - jkoQP - jkoP - kl - klQ 

+ j^oPQ + k^Q 



( 20 ) 



Since PQ = QP, we have the desired result. The first operator in (20) denotes 
an incoming wave (w.r.t. x) and the second an outgoing wave. We retain only the 
outgoing wave to obtain 

Pu = -ik^{Q - l)u (21) 

The square root operator Q is global in nature and we would like to make some 
approximations to derive a local operator from it. (It is global because when expanded 
in terms of series, it will contain terms of all derivatives). Now 



Q = 



2 I I da d I 

kl Q dy dy k^ dy^ ^ 



Let us rewrite Q as 



1 1/2 



Q = 



1 + [ii^(2/) - 1] + 

" ' 

<C1 normally 



1 da d 
a d{koy)d{koy) 



+ 



small 



d^ 

d{koyy 
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which may be expressed as 

c? = v/rrp 

w’here 

kla dy dy ^ kl dy'^ 
= a small operator! 



Treating \' as an algebraic factor < 1, we may derive the following rational approxi- 
mation {padc{] . 1)) 



Q = yr+T 



T — (Claerbout) 

4 + 3\' _ 2\' 

4 + r “ ^ ^ 4 + r 



so that 



Q-^ 



Substituting into (21) we arrive at 



21^ \ 
4 + \'J 



Pu = jkoiQ - l)i; 



Pxi = 



-2jkoV 
4 -f \ ' 



u 



or 



{4 + V)Pv = -2jkoVv 



. , 1 da 

using the fact that Q 2 — — we get 

Q dy 



/ 2 , - , d 






j ,,i a ■ 


(.1 + 3 )‘-o + <>=gJ + 5 ^ 


H 

II 


0 

O'! 

1 


{n + + 



The above equation is a wide-angle parabolic equation valid for propagation up to 
±20° [1]. Other approximations could be obtained by considering higher order pade 
approximants. 
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6. Boundary Condition on the Upper Boundary 



To truncate the computational domain, we consider a point high enough where the 
atmosphere is homogeneous with n = \ . The governing equation in the homogeneous 
region becomes 



Let us derive boundary conditions on a horizontal interface y = j/o- Tor the sake 
of simplicity we will derive boundary conditions of the mixed type on the interface 
y = 0 (instead of y = yo) given the initial data on the line x = 0 and boundary data 
on the line y = 0. Our derivation is based on the use of Fourier sine transforms as 
suggested in [2]. Although the basic philosophy of our approach coincides with that 
in [4], some of the details and the final results are slightly different from the latter. 



Consider the parabolic equation Uyy — 2jkuj. = 0 in a homogeneous region x > 0, 
y > 0, where k is a complex constant, 




subject to the initial condition u(0, y) = /(y). 0 < y < oo, and the boundary condition 
u(x,0) = y(x), 0 < X < oc. We assume that u(x,oo) — ► 0 and Vy{x,oc) — > 0. The 
equation 

Uyy = 2jkur (23) 

is of Schrodinger's type. We will treat the lossless case having a real value of k as 
the limiting case of the lossy problem having k ko — jc, e > 0. We will solve the 
problem using Fourier sine integral. 

Let 

[2 

f/j(x,A) = W-/ u(x,y)sinAydy (24) 

V 7T JO 
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Then using integration by parts, we see that 

^ Uyy{i,y)s\nXydy = ^ luy{x,y) s’m Xy 



y=0 



— Ai/(x, y) cos \y 



y=o 



u{x,y)s\nXydy 



Because i/y(j",cxD) — > 0 and u(x,oc) — > 0. we have 
I 2 yoo 



- Uyy{x,y) sill Xu dy = \ -Xu{x.O) - X‘^L\{x. X) 



TT Jo 



7T 



= \ -Xg{x)~ XH',{x.X) 



(25) 



Multiplying both sides of (23) with \J2/ ~ sin Ay. integrating over y = 0 to cc and 
making use of (25), we have 



^hg{x) - X^Usix, A) = 2jk^l\{x, X) 



or 



d 



A2 



A 2 



dx 2j k 

W’e may rewrite the above equation as 

dx 



= 5-r\ -»W 



2jl' V 






(2C) 



Replacing the dummy variable x in (26) with r and integrating both sides over r = 0 
to X. we arrive at 



b'.(7,A)e 



(A2/2j<:)7 



^ A 12 r 



~ / dx 

TT Jt=0 



or 



But 



t/,(x,A) = r 

Zjk \ TT 7r=0 

f7^(0. A) = u{0,y) sin Xy dy 

= \l~ /(y)sin Ayt/y = Fj(A) 
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A 2 r 



2jk V 7^ Jr=0 

Finally, taking the inverse sine transform on (27) we get 



u(x,y) = \ - [ s\n \yUs{x,\)d\ 

\ 7T Jo 



^(^)e(^V2ji)(r-x)^^ (27) 



= r F,(A)e-(^'/'^'>"sinAy(fA 

V 7T Ja=0 



+ — f A sin Ay / ^^drdX 

7T 2 7 A' 7a=o ^ Jt=0 ^ 



2jk 

2 



Defining 



= — [ [ /(r) sin At f/re sin Ay rfA 

7T J A=0 J r=0 

+ r r AsinAye('^/'^'=H^-")jA(fr 

Ti Ijk Jr=Q JX=0 

= — / /(t) / [cos(r — y)A — cos(r + y)A] (fr 

T Xt=0 ~'a=o 

+ - / vl / cos Aye^^"/^-’'‘'^(’'"^)(fA Jr 

7T Jt=o yA a ay j j\=o 

1 9 y • I \ 

= — f{x) I [cos(y — r)A — cos(y + r)A] ^^JA Jr 

7 T Jr=o Jx=o 

/\(x.y; xo.yo) = / cos(y — for ^ > ^o- 

Jo 



we write the expression for u{x^y) as 

1 



1 

u{x,y) = - /(r)[A'(x,y; 0,r) - A'(x,y; 0,-r)]cfr 

7T 7t=0 

We now evaluate the integral for A'. Consider 

l{a) = f cos , x > Xq 

Jx=o 



( 28 ) 
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Because of the exponential decay, we may differentiate under the integral sign to 
obtain 



— 1(a) = - r AsinQAc-<^'(^-^“)/2l'^l'X'-^'^“)JA 

da Jx=o 

j.-(A^(x-xo)/ 2 |fcp)(f-jfco) 



= r s\uQ. 

Jx=o 



d_ 

d\ 



(x-To)(c-jl-o)/|^'P 



d\ 






a|A-p 



(x-xo)(c-jAo)/|Ap 

r cos 
Jo 

jq|A-P r 

(j — To) A' 7o 

d/(a) ;qA 



=0 (^-^o)(( -jAo) 






(t - To) 



+ 



7(a) = 0 



da (x — To) 
or 

— = 0 /(a) = 7(q = o)c-^“'*^7(2(x-xo)) 

(7q ^ J 

Now 7(o = 0) can be written as 

/(a = 0) = r 

Let us view this integral in the complex A-plane. 



(29) 




Riit^'ow (S^ 
C^v^(^jla^ Ay<^(X) 




Rx(A) 



Complex A-Plane 
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For the integral to converge for {x — Xq) > 0, /?e[A^(e — > 0, i.e., 

> 0 =^> lm{XH'') > 0 



or 



or 



or 



Now 



0 < Arg(A2r) < 7T 
0 < 2 Arg (A) — Arg (k) < ir 



^ Arg {k) < Arg ^ ^ Arg {k) 



Arg {k) = - tan M -^ ) ~i- ■ -r < 1 



A’o 






~W, ^ ^ 1 - Wo 

I{q = 0) = / 

Jc 

where C lies in the region of convergence. 

In particular, let us choose a line from 0 to oc along the line C^: defined by 



A = 



I 2jk 

X - Xo 



in = 0 to oo 




4 



(30) 



(31) 



Zkr 



On this path, 



Arg(A)4 + lArg(.) = ^-^ 



l{a = 0 ) = 



2jk 



(^ - ^o) 



r 

J tiJ=0 



e-"” dw = 



wjk 



2(x - Xo) 



18 



^-jo'^k/(2(i-zo)) 



2{x - To) 



K{x,y; xo.yo) = 



TTjk 



2(t - To) 



^-J^(y-yo)^/(2(i-ro)) 



(32) 



Substituting this into (28). the field u{x,y) for the initial-boundary value problem is 
given by 



7:jA-Jt=o^ dy\\l2{x-T) j 

It is easy to see that 
0'‘K 



dy^ 



= r -X^ cos[(y - 
^0 



and, so for x ^ Xq 



d^< ^ OK _ dl< _ 1 d'^K 

dy^ ~ dx dx ~ 2jk dy^ 



which is the same equation satisfied by u. Now 

V In X J T=.0 



du 

a?*"-*'* 



y-0-^ 



Ijk 

27TT 






2jkr 



g-;/:rV(2x)j^ 






dr 



J/— 0+ 






5 



TTjk 



dr 



y-^0+ 



(33) 



(34) 
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in view of (34). Furthermore, we note that 



■^K{x,y; To, yo) = (joy; ^o2/o) 



I2jk 



Using this and evaluating the first integral by parts, we get 

' OO 

r=o 

+ - [ g{T)-^K{x,y: r,0) dr 

7 T Jr = 0 or ^ . 






y— 0+ 



TTX 

2 



ly— 0+ 



IX ^ 

i?(r)A'(j,2/;r,0) -f gr{T)K{x,0\ T.0)di 
Iv-o^ 

Jo 



2 

+ - 
7T 



= +J—fm + 

V TTX 



TTX 






l-rjk 2jk gr{r) 



2x 



7T 



r 9r{T) 

Jt= 0 yjx — T 



dr 



\‘ij± 

TTX 



[/(O)-y(O)] 



+1 






Jo y/x — T J 



7T [Jo \/x 

From the compatibility conditions on the initial and boundary values, we have 



Therefore 

du 



dy 



i^,y) 



limu(o',0) = 5(0) = lirnu(0,y) = /(O) 

X--+0 y—*0 



4= r r -^^dT 

y/x Jt=o Jt=o \/x — r 



y— 0 + 







( 35 ) 



It is easy to note from (35) that for k = ko — je, 



lim ^{x, y) 
T-^ody^ ' 



= 0 . 



y— 0+ 
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No matter how small the e is, we will use the same equality for the limiting case 
of c — ► 0. We will evaluate the integrals above approximately by replacing the 
derivatives with differences. 




We now consider the evaluation of dujdy\y=Q at x = Xp_i/ 2 - Consider initial data 
on the line x = 0. Let us assume that this initial data is known on a uniform grid 
t/„, = rnA?/. m = 0, 1,- • We approximate the derivative in the interval {ym-i- ym) 
by the forward difference formula 



dj{y) ^ Um - Um-i _ 
dy ^ Ay 



y ^ (j/m — 1 ) J/m ) 



where = u(0.r/m)- 



Then 



Let 



1 



•^"=0 m=l 

r"' 

y/x 



— 1 
Ay 
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We then have 



_L T'” = \!^ I 

yjx V K J_ 

y/kl(Tx)ym-\ 



yJkl(TTx)yn 






• y^*:/(rx)y m 
■ v/t/(7ri)y m — 1 



where F{fi) = complex Fresnel integral [8] 







(36) 



Consider now the boundary data on the tine y = 0. Assume that we have a non- 
uniform grid 0,Xi,X2) • • ■ )^p-i /2 = 2 ). On the interval (x^_i,Xto)> we approximate 
the derivative as 



9t{t) = 






1 

where u"’ = At the origin w’e have = uq. We write 

9t{t) ^ I^P-^ gr{T) I^P-U2 gr{T) 

— r Jt = 0 \/Xj,_i/2 — T Jr = Xp-i 



r 



\/-^P-l/2 



y/^p-ll2 ~ T 



:dT 



We evaluate the first integral as 



/. 



^P-i 5^(r) 



P“1 «,m «,m — 1 



\/a: — 



:c/r = 



^ r^rn 

Jxm-i 



= i; 



u’" - u’""’ 



— 1 



( — 2 y/x^Ti/T"— ~T ) 



zdr 



^m—1 






2 ^ ^ (>y^p-l/2 ^m-1 \/^p-l/2 

m = l ^m-1 



p-1 

2E 






\/5p— T/2 ^m — 1 ”1" \/^p— 1/2 ^ 



(37) 
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Now 



^P-)/2 gr{T) 



r 



y/Xp-\l 2 — T 



idr 



yV 1/2 _ yV 1 rlp-i /7 dr 



- u>'-‘ r^P- 

Xp_i Jxp^j 



{Xp 

= 2 



Xp~\f2 ^p—\ ] \J~X~p—\j2 "x 

yP-^!'^ - yP-1 /■Xp_i/2-Xp-l d\ 

-1/2 ~ ^p-1 ) ‘'0 \/A 



yP 1/2 _ j/P 1 



\JXp — \j2 ^p-1 

Substituting (3G). (37), and the above in (35), we get 



v/27 



du 


^P-U7 


dy 


y=0 


oo 




y] - Wm-l) 


m — 1 






- 



+ 



1 



8jt' 



x{Xp-^!2 - Xp-i) 



^p-\l7 

y=0 



\| 2^3:p_i/2 



ym - F 



\| 7TXp_i/2 



?/m-l 



^ m = l ‘^(^p — 1/2 ^m) “I" ^(Xp—i/2 ^m-l) \ 



+ 






^(^p— 1/2 ^p-l) 






(38) 



where 



D— 1/2 — r)(^p-l “I” ^p) ^ ^p-l/2 ^p-1 — r)(^p ^p-l) — r)*^^p— 1 



The above equation is the discrete version of a continuous boundary condition of the 
form 

N /X 

— + r( xjtv = s(x) 

dy 



where 



r(x) = 






Isjk 1 



y/ 3 : - 



p -1 



(39) 

(40) 








and 



_ 

^ fn=l \/^ “h \/x — 1 ^ 7 t(x Xp_j) 

F(x) = 

Jo 



u"-’ , (41) 



(42) 
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is the complex Fresnel integral [8]. 

For efficient implementation, the PDF and the various boundary conditions will be 
transformed into a curvilinear coordinate system generated by setting the lower ir- 
regular boundary as a constant curve curvilinear coordinate. 
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7. Transformations to a Curvilinear Coordinate System 



Consider the narrow angle parabolic equation with range dependent refractive index 
(oi 7 ^ 0) given in (18) 

Ur = 77-T; T S 0) + h r ^ 1 t)E 

together with the tropospheric boundary condition 

d 



Oy 



u{x^.yo) + r{xm)y{xm,yo) = ^(Xm.yo) 



an 



d the impedance boundary condition on the irregular boundary 

Ui/ C\U — 0 



■7 7 7 7 7 7 7- 









X-rY^j:)Jui. 

We will transform these to a curvilinear coordinate system 






c 



^ L- £ £ —C £ 




BC 



r 7- 





^ 








^ 1 






































' ; 




/ r 


/ ✓ 




^ 1 



c O 



Phvsical Domain 



Computational Domain 
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We assume that we have a transformation of the following form 

T = x(0, y = y{(,v) 

The various metrics needed in the transformed equations are [5] 



9u = ocl + y\ , gx 2 = + y^yn ■ 



We then have with ^ = x^yr, — Xj^y^ = 7^ 0 



Ux = ^(y^U£ -y^u^) = ^(u^y^ -u^y^) 

\/9 x^yr, 

1 / , X ■*^7? 

Uy = -—{-X^U^ + X^Ur,) = 

\/9 y^ 



^yy 



( 


1 


2/7777 




= — 


ci 

1 

i 


\yv / 


Vr , 





Substituting into the PDE we get 



1 



or 



Uf = 



^^yn 

x^a] 



{u^y^ - Ur,y^) = 



{2jk'o - a i) 



'^T] yriT) 

2 n 

yr, y-r. 



. I "77 aT]Tl 

+ a2- jUr, + 



^ {2jko-ai) 
Letting 



u + 



^2 3/777; 



X(_ 






.3/7; / ( 2 j^'o -Qi) yn. 



U-n + 



vl 



x^ 



{2jko - a i)y^ 



''T]7) 



b, = 



x^a\ 



{2jko — Gi ) 



b2 = 



fi 

yr, 



— 



O2 - 



X^ 



y^n 



+ 



y? 



y^ / (2;A’o - oi) x^ 



{2jko - ai)y2 



we express the PDE as 



Uj — 61U + b2‘Urj + b^U 



rjn 



(43) 



The normal derivative on a y = constant line is 

9\7 



Ui.(y = const.) = \ 7==^i 



9n 

^r? / 

9 vWn 






y( 



Xiyrj-yiXr, " ^xl +y|(x^y^ -y^Xr,) 



x^Xt, + y^y„ 



(44) 
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Variation of an arbitrarj’ vector r along the 77 = const, coordinate is 



= x^x + yi^y = 3(_ 

The unit vector along the tangent on y — const, is then 



s = 






^ y(.y 



l«d yjx\ + y\ + 



Defining 



cos 0 = 



sin 0 = 






yj 



we express (44) as 



y-y = 



{Xri cos <? + J/r 7 sin 6)\l(^ 



yr^cosO - Xr,s\nO 

COS 0 — Xr^ sin 9) 



With these substitutions, the boundary condition at the bottom boundary i 
0 gets transformed to 

(xr^cosd + yrjSmO) , 

Ur^ J ^ + [yr, cost/ — X Sill 0)CiU = 0 OU 7/ = C 

\l^l + yl 

For = 0 and using ^x^ + = x^/ cosO, we get 



V 

Uj, ^ sin cos + C)?/;, cos = 0 @ r; = 0 



Finally at the top boundary, we have 

^n/yv + ru = s 



or 

+ r{im, -^')yr,iim)H(m, N) = yrA(rn, ^’) 



1/ C]U — 

(45) 



(46) 



( 47 ) 
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8. Generation of Curvilinear Coordinate System 

Consider a piece-wise linear ground profile and a horizontal upper boundary 




Physical Domain 



Computational Domain 



^ — c — 









Z 

N 



^t+/ » ^1+1 



#-7 7 7 r - 



We generate the (x, y) coordinates of an interior point by using a linear interpolation. 
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Letting Ax, = (x,+, - x^), Ay, = (y,>i - y,), and A{, = ((f,+i - {,), we write 

Ax, 



X = X, + 



y 



A6 

n) 






I ^y^ [ t c \ 

y> + ~ ^') 



A6 



+ jyo 



(.<(< iM 

0 <T] < N 



( 48 ) 



At any interior point, (, < ^ < the various metrics are evaluated as 
Ax, 



A6 ’ 



x^ = 



y( = (i-jv 



X,, = 0 



(49) 



Ay, 

Af ’ 

^S1 



1 



Ay; / . • > 

yo-yt- ^(C - 6) 



yr,r, — 0 



At the boundary points, we use the central difference formulas to arrive at 

2'((( = <i) = 7 



- X, 


Ax,+i + Ax, 


(50) 




A^+i + A(f, 


-) 

n) 


Ay,+i + Ay, 
A{,+i + Ai^, 


(51) 



!/(({=«.) = 

Note that the analytical expressions yield a discontinuous value for these derivatives 
at the boundary points. We use the analytical expressions only to generate grid 
points and use the central difference formulas to arrive at the derivatives w.r.t. (. 
In other words, once the grid points are generated on the lines AA', BB', . . . , etc., 
we assume that the space is smoothly connected through the grid points. In the 
numerical implementation using Crank-Xicolson implicit scheme [6]. the metrics are 
needed at the midpoint w.r.t. (f. i.e., at { = {, + A^,/2 and the interior point formulas 
are applicable. For a uniform mesh in the computational domain. A^, = 1, y = y. 
y = 0.1,2,---,A’ 



X^ = 


Ax, 




yi = 


(-^) 


Ay, => 2/<(9 + 1 ) - 2 /«(?) 




1 ( 


2/1+1 + 2/1 ^ 


J/rj — 


N 


2 ) 




0 -^) 


( 2/1+1 + 2 /i\ , 9 


y = 


( 2 j + 


= 


2/1+1 + 2/1 
2 


+ 92 /r; = 2/0 + (9 — ■'^02/17 


) - y {(}) 


= yv 





Ay, 
' N 



(52) 

( 53 ) 

( 54 ) 



( 55 ) 
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9. Numerical Implementation by Crank-Nicolson Scheme 

Consider the narrow angle PE together with the boundary condition 

— b]U h2Ur) + b^Urjn (PDE) 

+ CiPr, COS = 0 at = 0 (BC2) 

Ur, + ryr,u = y„ s at y = N (BCi) 

We would like to implement the above using a Crank-Nicolson implicit scheme [6] 



sin 6 cos 6 






% 



t ' 



Ay\ 



- I 















P 



We use the notation = u{^p.rjg) = n{xp.yg) 

The various derivatives assuming = Ay — 1 are 

u^{ip^i/2,yg) = 

1 \ .1 



or 



Ur,[p--rqj = Up{^p.i/2,yg) ^ ~[u,^{ip_uVg) + UnUp,yg)] 



Ur,\p--,qj = - 



P— 1 p— 1 , p p 

- u,_, 



1 






( 1 A 


_ 1 


V~rV 


“ 2 



<-^1 +<+i - (<-l + <-i)] 
- 2uP + -f - 2ul~^ + 
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Substituting into the PDE, we get 









p-\ 



b2 



(<+} +<+i -<-i -^^^0 



^3 (p - 



Rearranging the terms we get 

1 /6 



(<+i - 2nP + - 2^^ ’ + 



bi' 



+ ^3 Ci-M+^3-7f < + 



2 V 2 



1 fb 






b^ 



<"^ + o(^3-^K:1 = o 



p-i _ 



Now let 






9 

P-1/2 






1 = 



p-l/2 



We then have for p = 1 , 2. • • •, q = 1 , 2. ■ • •. A' — 1 



QU 



9+1 



- (1 + 3)u^ + , = -quLJ - (1 - 3)u^„ ^ - 71 ^ 0-1 



(56) 



However, we will extend the applicability of this equation over q = 0, to 

accomodate the derivative boundary conditions on the lower and upper boundaries. 
For q = 0. we have 

anf - (1 + 3)ul + “ (1 “ ~ l^-V (5”) 

for q = A\ we have 

«^A’+i - (1 + + 7^^a'-i = -^>i^a"+i - (1 - (58) 
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From the BC 2 at the lower boundary we have 

V- 




— ^^sin0cos0^ (^0 ~ ^o~^) 



, , n^P-l/2 + ^0 \ n 

+ (ci2/„cos^?)P I ^ 1=0 

Multiplying this with and adding to (57) 



Letting 



(q + lOuJ - ^ - 27 

= -(q + 7)ur’ - - /3 + 27 

a = Q + 7 



Vv 

Cl j/^ cos 0 — 2 — sin 0 cos 6 



Uq 



C\y„ cos 0 + 2 — sin 0 cos 6 
^4 



u 



/3' = 0 - 27 ?/^ cos 6 ^ci - 
/3" = 13 - 2~,yr, cos 0 ( Cl + 



2 sin 
^4 ) 

2 sin 6\ 
^4 ) 



we may write the above equation as 

o'uf - (1 + /3')ul = - (1 - /?'>?■’ 

From the BCi at the top boundary we have 



< 

y 


> — i — i 

I 

1 


- 1 '- 




1 


^ 0 
» ^ - 


P-i 


^ 1 < 


u 

F 



' N'l 



(59) 
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073 



1 f p—1 p-1 p p ^ 

1 i ~ ^A'-l _i_ ^A' + l ^A^-1 ^ _|_ 



M 



16 jfc Uy + ^A’ ' 

7rAxp_i 2 



■J/'i 



— £( 
Ay 





/ 


k \ ( 


Wm-l) 




8ce 

3 

1 




. v\ 


^•^p-1/2 y 



it 



7TX„_ 



'Vm — 1 



P-1/2 



8A- g 



«A- 



+ 



7T 



= 1 1 /2 •^m \/-^p — 1 /2 ^tt4 — 1 ) \ 



16jl- 

7rAXp_, 



,p-i 



or 



»^A’+l + ^^A- + l - (^A-'l + »'A-l) + ypS 
multiplying with — Q^- and adding to (58) 

1 + /? + yp8o 



jA’ 



7T^ Xp_ ] 



(ufv + u^v ’) = yp45P 



y'A- 



\ 7rAXp_i 



l^A- + (*■ + o)u^v_l 



1 - /5 - yp8o. 






\ ttAx 



p-1 



^^A- ’-(')+ o)u^v-] - 4a?/p5P 



Letting 



\ = 3 -\- ypSo 






/. Axp_i 



')' = O+O 

we rewrite the above equation as 

(1 + A)ufv - l'ufv_i = (1 - A)u5,“' + -)'ufv_', + iay,s’"'l'‘ 



(60) 



where 

cP-l/2 



Ay h 





/ 


k \ ( 


^m- 1 ) 




1 

6 

Si 




. U 


7TXp_i/2 y 



\ 7!'Xp_i/2 



ym -1 



-)?£ 



7,™ _ 7/'"“J 

UA' U;,,. 



^ m=] ^\/2'p— 1/2 "t" ^Xp_i /2 2^m— 1 ^ 



+ , 



16yl- p_, 

7rAxp_, 
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and 



1 



a^p-i/2 



Xp Xp^\ 

1a 1/ 

^p— 1 d” 2 ' p ^p~i) 



We augment the equations given in (56) for g = 1 , 2, . . . , A' — 1 with (59) and (60) for 
7 = 0 and A’, respectively to define for 7 = 0, 1,2, . . . , A^ The system of equations so 
defined can be expressed as a matrix equation of the form 



■ 


A 






’ ^0 ' 




■ A^ 


A 


- 




r 1 




A^ 


A 








A 


A 


A 








A" 


A A 










A 


A A 










A A _ 




. - . 








1 




. . 



0 



+ 



( 61 ) 



iay^s” _ 

where A' denotes a non-zero entr\’. The tridiagonal matrix on the left hand side of 
(61) can be inverted efficiently to yield a solution on line ^ = i^p in terms of the field 
values on the line ^ — (p-\. Equation (61) can be used to march forward in range 
starting from initial data specified on if = 0. 
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